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ABSTRACT 



The small perturbation, two-dimensional transonic equation 
is manipulated with a separation- of -variables approach to 
obtain two ordinary, nonlinear, differential equations. 
Numerical integration of these differential equations results 
in new transonic boundary surfaces for planar external flows. 
A key ingredient in these solutions is the identification of 
dependence of two integrations constants, a and /?, on the 
parameter (1-M^ 2 ) . The anticipated behavior for both the Mach 
number and the pressure coefficient is used as a guide in the 
actual selection of the adjustable constants in the problem. 
The physical reality of our boundary surfaces is examined by 
displaying the boundary conditions they satisfy. The strictly 
sonic flow (M oo =1.0) has an analytic representation 
corresponding to a divergent surface which goes supersonic. 
This sonic solution is compared with an Euler-CFD approach 
confirming the validity of our results over the region where 
small perturbations apply. Solutions are also shown for 
M„=0 . 8 , 0 . 9 , 1 . 1 , and 1.2 . These results are consistent with 
known behavior for both subsonic and supersonic external flow. 
Since the results of this work yield actual transonic 
contours, we can examine shockless surfaces for design 
applications. The possibility of starting with transonic 
surface is of interest to present day CFD approach. Finally an 
entire transonic upper surface is presented for M oc =0.8, by 
patching a subsonic Mach number, which reaches a plateau at 
M=1.0, with a sonic flow. This patching requires the careful 
interpretation of a nondimensional reference length, called 
Y 0 , which is a function of M w . 
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I . INTRODUCTION 

The governing nonlinear equations, together with realistic 
boundary conditions have always been the main difficulty in 
transonic aerodynamic analysis, [Ref. 1] . While the hodograph 
transformations approach linearizes these equations, the 
transformations are confined to a limited number of transonic 
flows due to complicated boundary conditions. The small 
perturbation approach simplifies the full potential 
Equations; furthermore, like other predictive methods, it 
succeeds in retaining the great sensitivity of transonic flows 
within small perturbations. This paper expands on the exact 
solution of the small perturbation, nonlinear, two 
dimensional, transonic equation, using the guidelines and 
tools given by Biblarz 0., [Refs. 2 and 3]. Starting with a 
separation- of -variables approach, Chap. Ill reveals two 
nonlinear, ordinary, differential equations which lead into 
two exact solutions. Chaps. IV and V describe the numerical 
integration of these implicit solutions, which finally yield 
the boundary surfaces. Later, these boundary surfaces are 
found to satisfy the boundary conditions for a two dimensional 
surfaces. These findings are then compared with a finite 
difference solution of the Euler Equations for the sonic 
solution. The usefulness of the shockless, transonic, boundary 
surfaces appears once transformed to dimensional body 
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surfaces; "patching" (i.e., translating), the dimensional, 
subsonic, body surface of C p =0 at the upstream inflow, with 
the dimensional, sonic body surface, enables us to form 
complete transonic upper surfaces of interest for design. 
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II. 



THE TRANSONIC EQUATION 



The small perturbation nonlinear, two dimensional 
Transonic Equation 

, ^ M 2 (y+1) ^ ^ ... 

(1-M. 2 )^ + <byy = &X&XX (1) 



is described by the velocity potentials 

u = C> x v = $ y (2) 

Multiplying eqn. 1 by 1/U^ enables us to transform the 
Transonic Equation to 

(I-mZ)^ + = M^y+D^xx < 3 > 

where the nondimensional velocity potentials are 



4>x = 



u 

U m 






v 

~U m 



(4) 



The Meyer solution for the de-Laval nozzle defines a 
velocity potential 



4 > (x,y) 



= b 



— + b(y+l)Mi^ +b 2 (y + 1) 2 Mi^~ 
2 2 24 



(5) 



which satisfies the nondimensional form of equation (3) for 
M oe =1.0 . The Meyer solution belongs to a general class of 
sonic flows as reported by Guderley and Yoshihara, [Ref. 4] . 
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The transonic equation (1) may be modified to a more 
useful form if multiplied by (y+1) /U„ , instead of just 1/U- 
«, to give 

+ 4> yy = (J)^ (6) 

Here the modified velocity potentials are 

4>* = Mi (y + 1) -jj- 4> y = M„ 2 (y + 1) (7) 

For a given problem, and y are constants, and eqn. 6 is 

really less cluttered. However, notice that M„ J (y+1) /U„ is 
already given in the derivatives of 4> in eqns . 7, and 
therefore appears in the expressions of C p and the boundary 
surfaces as will be shown later. This form of the modified 
transonic equation (6) will be examined in the next chapters. 
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III. 



THE GOVERNING DIFFERENTIAL EQUATIONS 



An exact solution to the modified transonic equation 6 has 
been given by Biblarz, [Ref. 2], starting from 



where $ s (x,y) refers to the sonic, small perturbation solution 
(M w =1.0) . In other words, "any flow which satisfies the sonic 
equation also satisfies the full transonic range with similar 
boundary conditions," [Ref. 2]. 

An exact solution may be obtained by using the separation- 
of -variables approach with the potential function <£(x, y) of 
the form 



Substituting the above 4 > function in the modified 
transonic equation 6 results in two ordinary, second order, 
nonlinear differential equations 



4 >(x,y) = 4> s (x,y) + ( 1 -Ml ) x 



(8) 



<J >(x,y) = £ (x) r| (y) + (1 -Ml) x 



(9) 




(10) 




( 11 ) 



where X is a separation constant. 
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A simple solution to the first differential eqn. 
obtained by multiplying both sides by d£/dx, 



dx\dx dx 2 ) 



- (AO = 0 

dx 



or 



- Af 2 

c?x[ 3\ dx) 2 ^ 



0 



Thus 



dt = 3 1 lA f2 

dx N 2 



+ a 



dx = 



dt 



\ 



3 A z2 
2 



+ a 



and 



x - x n 



-!■ 



dt, 






3A 



o + a 



where a and x 0 are integration constants. 

A similar procedure can be performed to solve the 
differential eqn. 11 yielding 



cfrf d M 

dy( dy 2 ) 



PP (Ari 2 I 
dy 



= 0 



or 



d F 1 / dr| \ 2 
dy[ 2\ dy) 




= 0 



10 is 



( 12 ) 



(13) 



(14) 



(15) 



(16) 



second 



(17) 



(18) 
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^P- = ±\ — n 3 + P 

dy N 3 1 K 



( 19 ) 



dy = t gn 




( 20 ) 



and 



y ~ Yo = ±/-=^ 

2A„ 3 



( 21 ) 



+ P 



where (3 and y 0 are integration constants. 

Eqns . 16 and 21 are exact solutions to the differential 
eqns. 10 and 11. In fact, boundary conditions are needed to 
evaluate the integration constants a, |3, y 0 as well as X. 
An important step in this procedure is the recognition of cn 
and /3 as functions of (l-M,,, 2 ) , [Ref. 2], i.e.. 



The upper positive sign and the lower negative sign in the 
above equations will be shown to belong to M^al.O and M oc <1.0 
respectively, namely, 



a = ±q (1 -Mi) 2 



( 22 ) 



p = ±C 2 (1 -Mi) 



( 23 ) 



aaO & 0 for 



Mo.al.0 



a<0 & /?< 0 for 



M OJ <l . 0 



See Appendix B . 
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The constants, a, j3 and X can be shown to be related by the 



expression 




( 24 ) 



Since , C 2 and X are positive constants 



C,C, = — X 
1 2 3 



( 25 ) 



The usefulness of these definitions will be outlined in 
more detail once we use the exact solutions from eqns . 16 and 
21, while verifying the correct boundary conditions for the 
pressure coefficient C p . 

Upon inserting eqns. 22 and 23 in eqns. 16 and 21 we have 




( 26 ) 




( 27 ) 



Factoring out C, and C 2 we obtain 




( 28 ) 
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_ _1 

y - y 0 = ± c 2 2 



/■ 



dr 1 



N 



2X 
3 C, 



t) 3 ± (l-wi 



Introducing new dependent variables 



I 



3 A ^ 



2 C. 



i/ 



*1 



= T1 




1 

3 



the above simply becomes 



x " * 0 





dl 

± ( 1 -Ml) 2 



y - y 0 




/ 



^rj 3 ± (1-Mi) 



Finally, if we define new independent variables 



X = {x 





1 

2 



Y = 



i 

(y - y 0 ) c 2 6 




(29) 



(30) 



(31) 



(32) 



(33) 



(34) 



(35) 
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eqns. 32 and 33 then follow as 



X = 




_dl 

± (1 -Ml) 2 



( 36 ) 



Y = 




dT) 

\Jf\ z ± (i-mI) 



( 37 ) 



This set will be numerically integrated to obtain 
transonic surfaces with their corresponding C p and Mach number 
profiles . 
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IV. 



NUMERICAL INTEGRATION 



We can rewrite eqns . 36 and 37 in a more general form and 
specify the two integral limits 

±1 



As indicated before, the upper and lower sign within the 
integrand root and the integral limits is used for M^.al.0 and 
M co <1.0 respectively. In addition, our flow of interest is only 
the first quadrant in X and Y, which represents physical 
external flows. To show this clearly, for the subsonic case 
M„ < 1 . 0 . eqns. 38 and 39 reduce to the form 




o (z 2 ± (1-Aff) 2 



( 38 ) 




( 39 ) 







1 
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and for M„al.O we have 



X=) - UO ( 42 ) 

o sjz 2 + {1-Mi) 2 



Y = 



- / 

2 

-d-«f) 3 



dW 

(w 2 + (1-Afi) 



J. 

rj^-(l-Ad) 3 



( 43 ) 



However, we will frequently retain the more general "±" 
notation in our discussions further on. 

Before we head to numerical solutions, we should indicate 
that sonic flow (M co = 1.0) can be given in explicit form through 
eqns. 40,41,42,43, i.e., the resulting equations for sonic 

flow are 

X = f ( 44 ) 

J o 'jz* 



Y = 



f JM. 

o 



( 45 ) 



Solving equations 44 and 45 we obtain our sonic solution 




( 46 ) 



= 




( 47 ) 
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Having defined %,rj,X and Y in eqns . 30,31 and eqns . 34,35, 
we can express the sonic solution rewriting eqns. 46 and 47 as 




( 48 ) 



( 49 ) 



Providing the sonic perturbation potential from eqn. 9 



4> s = 5r, 



( 50 ) 



and inserting eqns. 48 and 49 we have, [Ref. 2] 



4> s 



1 U-x 0 ) 3 

3 (y-y 0 ) 2 



( 51 ) 



Next, eqns. 40 to 43 are numerically integrated and 
plotted in Figs. 1 and 2 for various , (e.g., 
0 . 8 , 0 . 9 , 1 . 0 , 1 . 1 , 1 . 2 ) . Plotting rj vs. Y for Y>0 in Fig. 2 is an 
appropriate choice for investigating physical external flows 
as stated earlier. 

Our numerical solution has been based upon the Newton- 
Cotes method of order 4, [Ref. 5 ] , which may be applied to an 
even number of subintervals, each described as, 

b 



ff(x) dx 

a 

( 52 ) 



(b-a) 

90 



7 f(a) +32 




(a+b) 

4 



+7 fib) 



(See appendix A) . 
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Numerical solutions which express | and r? , as has been 
noted, can be combined with eqns . 30,31 and eqns . 34,35 and 
substituted in eqn. 9 to give a complete solution for the 
perturbation potential in terms of the positive constants X or 
(C, , C 2 ) , and (xq, y 0 ) . 
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V. 



BOUNDARY SURFACES 



For an inviscid flow, the condition to be applied at the 
surface of a solid boundary is that the direction of the flow 
velocity vector is tangent to the solid surface, [Ref. 6] . In 
terms of perturbation velocities this boundary condition 
becomes 



Idy 

\ dx 



surface 



V 

CL 



(53) 



The modified velocity 
earlier is 



perturbation potential as pointed out 



4> y = Afi ( y + 1 ) -jj- 

^ oo 



(54) 



thus 



v _ 4y 

u - Mhy+D 



(55) 



We can show by substitution that 

l_dy\ = 4 V 

l I surface M? ( Y+l) 



(56) 



Now, deriving the perturbed potential function in eqn. 9 
with respect to y gives 




( 57 ) 
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Eqn. 19 for y>0 becomes 



2A, + p (58) 

dy 1 3 1 H 



This relation can be rearranged using rj and /3 in eqns . 
23 to yield 



= -< 

dy 



C 2 2 sj f) 3 ± (1-Mi) 



31 and 

( 59 ) 



Upon using this equation and the relation of eqn. 30 to the 
derived velocity potential in eqn. 57 results in 



4> y = -l 




d-Mi) 



( 60 ) 



The last equation can be rewritten with the benefit of the 
related constants C, , C 2 and A in eqn. 25, hence 



4> y = -|?7 rj 3 ± (1-Mi) 



( 61 ) 



Once we have found the form of </> y , the expression for the 
surface boundary condition is given as 



11 


_ 2 


1 


u 


\ dx 1 surface 


3 


Mi (Y+l) 


^ V 



fj 3 ± (1-Mi) 



(62) 



However, an important result to observe is that | and r) in the 
R.H.S. of this equation are given implicitly as a function of 
X and Y respectively in our earlier numerical integration 
solution, while the L.H.S. is expressed by x and y. To get 
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compatibility, we then transform the L.H.S. of eqn. 62 using 



the chain rule 



dy _ dY dX/dx 
dx dX dY/dy 



( 63 ) 



Deriving eqn. 34 with respect to x and eqn. 35 with respect to 
y gives 



dX = "•?/ 3A\-| 

dx 1 \ 2 j 



( 64 ) 



dY _ r T/2 A.\t 
dy M 3 J 



( 65 ) 



therefore 



dy = 3 dY 
dx 2 dX 



( 66 ) 



Equation 62 can be rewritten using the above transformation as 



dY 

dX 




1 

Ml ( Y + 1) 



l)/ Tj 3 ± (1 -Ml) 



( 67 ) 



or 



dY 


= - ( 2 \ 2 1 


J rj 3 ± (1 -Ml) 


' 3 / mI ( y + i ) 



l dX 



( 68 ) 



We are now in a position to compute the exact 
surfaces out of eqn. 68. Starting with the L.H.S. 
define 

E(Y) = 1 

fj 3 ± (1-AL 2 ) 



boundary 
we then 

( 69 ) 
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Having defined eqn. 69 we can plot E(Y) vs. Y in Fig. 3 for 
various M„, (e.g., 0.8,0. 9,1. 0,1. 1,1. 2), keeping in mind our 

"±" notation for the supersonic (upper sign) and subsonic 
(lower sign) . 

The asymptotic value of E(Y) occurs at Y 0 where 

fj - * (1 -Ml)' 1 (70) 

In other words, the asymptotic value can be expressed as 
Y 0 =Y 0 (M W ), or as function of (1-M„ 2 ), which has been 
numerically found to be 

Y - 2 -4 

0 i (71) 

|i-Mi | 6 

Y 0 vs. is sketched in Fig. 4 to demonstrate an interesting 
symmetry with respect to M oc = 1.0 . Here Y 0 =oo but note that for 
practical reasons we can adopt finite values of Y 0 (perhaps as 
low as 7 to 20) to represent conditions sufficiently close to 
M„=1.0 . 

Defining the integral of the R.H.S. of eqn. 68 as 

X 

K(X) = -(4) 2 __! f l dx (72) 

\ 3 / Y + l) o 

and redoing numerical integration allow us to plot K(X) vs. X 
for 7=1.4 in Fig. 5. It is important to realize that K(X)>0 
for M w <1.0 and K(X)<0 for M„al.O , before solving eqn. 68. 
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Therefore, integrating both sides of eqn. 68 yields for K(X)>0 

i'o 

f E(Y) dY = K(X) (7 3) 

Y 

and for K(X) <0 

y 

(e(Y) dY = K(X) (74) 

i'o 

Once we specify X for K(X) we can immediately solve the L.H.S. 
integrals of eqns . 73 and 74 for Y, expressing 

(Dsurfac* = YW ( 75 ) 



namely, the required boundary surface for any chosen . 

We can rewrite eqns. 73 and 74 in the more complete form 
for M 0o <l . 0 as 



I'd 



/ 



dY 

yj f) 3 - (1-Ad) 



(If 

\ 3 I Mi ( Y + 1 ) 



h ^ 



0 



and for M„£l . 0 



/ 



^0 



dY 

/ n 3 + (1-Mi) 



Ilf 

' 3 ' Mi(y+1) 




(76) 



(77) 



Eqns. 76 and 77 are equivalent to eqns. 73 and 74 which 
finally yield the boundary surfaces. These are shown 
nondimensionaly in Figs. 6 and 7 for M oo = 0 . 8 , 0 . 9 , 0 . 97 and in 
Fig. 8 for M„=l . 0 , 1 . 1 , 1 . 2 . 
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The sonic surface is most easily verified by direct 
analytical expansion of eqn. 77, using the sonic relations 



stated in eqns . 46 and 47. 




-3 



dY = - 



2 \ 2 1 
3/ (y+1) 




x 



( 78 ) 



Solving the above integrals on both sides gives 



y 4 - y 0 a = - 




( 79 ) 



Recall that Y 0 -» oo for the sonic case. However, choosing for 
example Y 0 =20 gives an asymptotic value of ij=0.01, (eqn. 47), 
and E(Y 0 )=1000, (eqn. 69), for our explicit sonic solution. 
(See also Figs. 2 and 3). Therefore, dividing by a finite Y 0 
enables us to obtain the nondimensional sonic boundary surface 
as a function of y • 



This will be shown to be a divergent surface which increases 
the Mach number to supersonic conditions. 




( 80 ) 
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VI. 



PRESSURE COEFFICIENT 



For a flow with small perturbations the 
coefficient is given by, [Ref. 8] 



— + (1-W„ 2 ) 
U 




A first order approximation for the linearized 
coefficient in a two dimensional planar flow is 

r = - 2u 
p U 

oo 

Recall the modified velocity potential 

<|> x = Y + l)-^ 

^ CO 



Thus 



u_ = 

Mh y+1) 



and from eqn. 82 we see that 



C p = - 



2( t>> 



M 2 (y+1) 



Now, deriving the potential function in eqn. 
respect to x we have 

4> x = + (1 -Ml) 

x dx 



pressure 

( 81 ) 

pressure 

( 82 ) 

( 83 ) 

( 84 ) 

( 85 ) 

9 with 

( 86 ) 
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Rewriting eqn. 14 using the relation for a in eqn. 22 as 



= 

dx \ 



^Z 2 ± C 1 ( 1 -M 2 ) 2 



( 87 ) 



By taking advantage of £ as stated in eqn. 31 and factoring 
out C 1( eqn. 87 becomes 



-^ = q 3 7 £ 2 ± (1-Mi) 



( 88 ) 



Substituting eqn. 88 and the relation for rj in eqn. 31 to the 
velocity potential in eqn. 86 yields 






( 89 ) 



Once again, using the benefit of the related constants C, , C 2 
and A in eqn. 25, eqn. 89 becomes 



<\> x = rj 7 V ± (1-M^) 2 + ( 1 -A/f ) 



( 90 ) 



therefore eqn. 85 for C p follows as 



CL = 



-2 



MZ (y + D 



i) 2 {V~± (1-M^) 2 + (1-M» 2 ) 



( 91 ) 



We must anticipate the above expression for C r to satisfy the 
boundary conditions on external boundary surfaces as will be 
explored further on. 
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VII. 



THE LOCAL MACH NUMBER 



The velocity vector V is given by the undisturbed uniform 
velocity and the perturbation velocities for the 2D flow as 

V = £(U„+u) + jv (92) 



Thus, we may obtain the local Mach number M in terms of the 
above velocities and the local speed of sound "a". 



M z = 




(U„ + u) 2 + v 2 



(93) 



Neglecting higher powers of the perturbation velocities gives 



M 2 



(Cl 




a 



2 



(94) 



In addition, the ratio of the local speed of sound to the 
undisturbed uniform speed of sound becomes 



a 2 



a 



2 

00 




(95) 



or 



a 2 = T 0 /T„ 
2 T / T 

<3co i 0 / 1 



(96) 



were T 0 designates the stagnation temperature. 
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Eqn. 96 can be replaced using the well known relation 



1 + X—LmZ 

2 

1 + XZ±m 2 
2 



( 97 ) 



Substituting the local speed of sound from eqn. 97 to eqn. 94 
results as 



M 2 = AC 



1 + 



2 u 
U„ 



‘ Y-l 2^ 
1 + -±— -Ml 
2 

1 + XzX-M 2 

2 



( 98 ) 



Having stated C in eqn. 82, eqn. 98 is reduced to the form 



M 2 



i + 1—Lm 2 

2 



AC (1-C p ) 

i +1z^mI 



( 99 ) 



Rearranging eqn. 99 and factoring out M : yields the final 
expression for the local Mach number 



M 2 



AC (1-C p ) 



1 + X — - 

2 



mIc„ 



( 100 ) 



Note that M=M 0C when C p =0 as it should. 

This important relation, combined with the corresponding 
eqn. 91 for C p , will be used when discussing the boundary 
conditions in the next chapter. 
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VIII. 



BOUNDARY CONDITIONS 



Rather than providing the boundary conditions for the 
governing differential eqns . 10 and 11 to evaluate the 

constants, we imposed a solution by defining a and /3 in eqns. 
22 and 23 respectively. That is to say, boundary conditions 
have been already implied in the definitions of a and /3. This 
procedure led us to reveal the boundary surfaces upon using 
the tangency condition and compute the pressure coefficient C p 
on the body surface, (eqn. 91) . 

Thus, we arrive at the necessity to verify the validity of 
the expression for C p for the boundary surface, and verify the 
following requirements: 

• At the upstream end (X/Y 0 =0) , we must require that C p =0 
ensuring the unperturbed Mach number of the transonic 
inflow to equal M B . 

• For the subsonic inflow at the downstream end (X @ Y/Y 0 =l) 
we must require that the local Mach number saturates 
towards M=1.0 . This is because of the well known relation 
in gas dynamics for a convergent flow. 

We first restrict our attention to the first requirement, 
for M^al.O at X/Y 0 =0 

l = = o (101) 



At Y/Y 0 =l eqn. 70 becomes 

i 

n - n* = -d -mI) 3 



( 102 ) 
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Inserting the above terms of £ and rj in eqns . 90 and 91 gives 



and 



4> x = -(1 -Ml) 3 7 (1-Mi) 2 + (1-wi) = 0 



( 103 ) 



C @ ( X/Y 0 = 0 , Y/Y 0 = 1) = 0 



( 104 ) 



Therefore, from eqn. 100 

M® ( X/Y 0 = 0 , Y/Y 0 = 1) = ML (105) 

Hence, eqns. 104 and 105 comply with the first requirement. 
The complete behavior of the pressure coefficient and the 
local Mach number are plotted in Figs. 9 and 10 for M oo =1.0 and 
M oe = l.l. In addition, it is important to note that the boundary 
surfaces ought to be truncated at M^sl.2, so as not to exceed 
the small perturbation assumptions. 

Now, for M„ <1.0 at X/Y 0 =0 

l = l* = 0 (iOS) 

However, at Y/Y 0 =0 

f=j - fj ' = oo ( 107 ) 

Thus, eqn. 91 results as 

C p @(X/Y 0 = 0 , Y/Y 0 = 0) - ~ (108) 

We immediately conclude that the above coordinate (0,0) can 
not serve as the upstream start for the subsonic inflow. 
Instead, we must localize a new upstream end coordinate along 
the boundary surface whose £* and i)* values satisfy eqn. 91 
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for the equality C p =0 . Since every coordinate on the boundary 
surface (X/Y 0 ,Y/Y 0 ) is a function of | and rj , (eqn. 76) , we can 
numerically scan the locus of all the coordinates, i.e., the 
boundary surface, for such values of f* and rj * that will 
satisfy eqn. 91. For example, scanning the results for the 
subsonic boundary surface of M oo =0 . 8 reveals; £ = £*=- 0 . 35763 and 
rj = rj*~3 .01393 at the coordinate (X/Y 0 =0 . 313 66 , Y/Y 0 =0 . 40616 ) . 
Inserting the above values of £* and rj* in eqn. 91 yields 
C p =0.0003 . Hence, by relocating a new upstream coordinate for 
the subsonic inflow we are then able to satisfy our first 
requirement. (This is shown later in Figs. 11 and 12). 

We now consider the second requirement for the subsonic 
outflow at the downstream end where 

l - V = -(1- Ml) 

1 

rj @(r/r 0 = i) - fj* = d-Mi) 3 

or equivalently dY/dX -* 0 
Substituting these terms for £ and rj in eqns . 90 and 91 gives 

4> x = (1 -Mi) (111) 



(109) 

( 110 ) 



and 



-2 (1-Mj) 
Mi (y+1) 



( 112 ) 
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Thus, we can express the local Mach number using eqn. 100 




Ml (y+1) 



( 113 ) 



Upon expanding eqn. 113 and factoring out M„ 2 we obtain 



Hence, we have verified the second requirement: The local 
Mach number saturates at the downstream end to the value M=1.0 
for any subsonic inflow M oo <1.0 . The pressure coefficient and 
the local Mach number behavior along the boundary surface are 
plotted in Figs. 11 and 12 for M„=0.8 and M„=0.9 respectively, 
showing clearly both the upstream and downstream boundary 
conditions . 

Knowing that boundary conditions require the gradient of 
the potential function 4> to vanish far away from the body 
gives rise to an important question: How is the flow pattern 
along the boundary surface propagated into the flow field? We 
treat this question by considering a computationally efficient 
flow solver developed by J. A. Ekaterinaris , [Ref. 7], This 
CFD solver consists of the inviscid Euler equations and uses 
fully implicit two dimensional Crank-Nicholson method with a 



Ml (y-1) +2 
Ml (y-1) +2 



( 114 ) 



or equivalently 



M = 1 



( 115 ) 
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central difference ADI scheme, and fourth and second order 
numerical dissipation, [Ref. 9] . This solver has been modified 
here to solve the inviscid flow around our earlier boundary- 
surface solution for the chosen sonic case M„=1.0 . By- 
creating a grid mesh which consists of a boundary surface 
edge, we are able to map the flow pattern. Upon observing the 
constant Mach lines as graphed in Fig. 14, we deduce that the 
flow is shockless through the whole transonic range; a shock 
is formed further downstream at the local Mach number M=1.9 . 
These results are in complete accordance with the inviscid 
small perturbation transonic solution, (see Fig. 9) . In 
addition, each vector on Fig. 13 describes the resultant flow 
direction at a specified grid point , (not the streamlines !), 
colored by Mach levels. Thereby, we may evidence the far- field 
boundary conditions 

V=C7„ u = v = 0 (116) 

In other words, the flow pattern is propagated out from the 
boundary surface gradually towards satisfying the far- field 
boundary conditions far away from the body. 

It is of special interest to inspect the dimensional 
boundary surface for possible design applications; For 
example, the dimensional, subsonic, boundary surface of 
M oo =0.8 is obtained by multiplying the nondimensional 
coordinates (X/Y 0(08) , Y/Y 0(08) ) by Y 0f08) =2 . 84 , i.e., the 
dimensional, downstream end coordinate results as (X, ,2.84), 
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15a. 



as plotted in Fig. 15a. Similarly, multiplying the 
nondimensional , sonic coordinates (X/Y 0(1 0) , Y/Y 0(1 0) ) by Y 0(10) =14, 
yields the dimensional, sonic, boundary surface, where the 
dimensional, upstream, starting coordinate becomes (0,14). 
Now, translating the dimensional, sonic, boundary surface, 
(Fig. 15b), by (AX, AY) = (X, , - ( 14 - 2 . 84 ) ) , to match the 
downstream end coordinate of the above dimensional, subsonic, 
boundary surface, "patches" both dimensional boundary 
surfaces, to form a 2D planar upper surface as shown in Fig. 
15. We may use this technique to "patch" any subsonic boundary 
surface in order to form our transonic upper surface of 
interest . 
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IX. 



CONCLUSIONS 



We have shown that an exact solution for the transonic 
equation satisfies the complete behavior of the pressure 
coefficient, and the local Mach along a 2D- surface for planar, 
external flow. It is of special interest to inspect the 
dimensional boundary surface for possible design applications; 
For example, the dimensional, subsonic, and sonic boundary 
surfaces are obtained by multiplying the nondimensional 
coordinates (X/Y 0 ,Y/Y 0 ), by the corresponding Y 0 . Translation 
of the dimensional, sonic, boundary surface, to match the 
downstream end coordinate of the above dimensional, subsonic, 
boundary surface, "patches" both dimensional boundary 
surfaces, to form a 2D planar upper surface as shown in Fig. 
15. This technique may serve as a powerful tool for the 
implementation of supercritical, shockless, planar surface 
designs, because of the saturation of rj with Y (Fig. 2), we 
can also describe certain internal flows such as planar 
converging, diverging tunnel walls. We may also extend these 
results to axisymetric transonic flow. 
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APPENDIX A 
COMPUTER PROGRAMS 



* THIS FORTRAN PROGRAM USES NUMERICAL INTEGRATION BASED ON 
THE NEWTON-COTES FORMULA TO SOLVE NUMERICALLY EQNS . 40 AND 42. 

THE OUTPUT IS IN THE FORM OF ^ vs. X 

************************************************************ 

program KSI 
PARAMETER (N=5) 

IMPLICIT REAL *8 (A-H,0-Z) 

IMPLICIT REAL *8 (M) 

IMPLICIT REAL* 8 (K) 

INTEGER INDEX, I 

DIMENSION X (N) ,XDOT(N) ,SAVEX(N) , SAVED (N) ,KSI(N) 

DIMENSION MACH (N) , XI (N) , X2 (N) 



MACH ( 1 ) =0 . 8D0 
MACH (2) =0 .9D0 
MACH ( 3 ) =1 . 0D0 
MACH ( 4 ) = 1 . ID 0 
MACH (5) =1.2D0 



KSI (1) =0 . 0D0 
KSI (2) =0.0D0 
KSI (3) -1.0D-3 
KSI (4 ) =0 . 0D0 
KSI (5) =0.0D0 

C MAIN PROGRAM * 

c*********************************************************** 

DO 100 1=1,5 
IF ( I . LT . 3 ) THEN 
H= - 0 . 001D0 
GOTO 5 
END IF 



5 
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50 



H=0 . 01D0 

IF ( I . LT . 3 ) THEN 

TF= - ( 1 . OD 0 - MACH ( I ) **2) 

GOTO 10 
END IF 
TF=12 0 . 0D0 

I F ( I . LT . 3 . AND . KS 1(1) .GT.TF.OR. 

I.GE.3 .AND. KSI (I) .LE.TF) THEN 
do 50 index=l,5 

CALL DERIVative (KSI , N, I , XI , X2 , XDOT, MACH) 
CALL NEWTONCOTES (KSI , XDOT, X, SAVEX, SAVED, H 
end do 



N, INDEX, I) 
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WRITE (10+i, *) X ( I ) , KSI ( I ) 

GOTO 10 
END IF 
100 END DO 
STOP 
END 

c* ************* * ********** * * * * * * * ********************* * * * * * + 

c SUBROUTINE DERIVative * 

£*********************************************************** 
SUBROUTINE DERIVative (KSI , N, I , XI , X2 , XDOT, MACH) 

IMPLICIT REAL* 8 (A- H, 0- Z) 

IMPLICIT REAL* 8 (M) 

IMPLICIT REAL* 8 (K) 

DIMENSION XDOT (N) , XI (N) , X2 (N) , KSI (N) , MACH (N) 

GO TO (20,20,10,10,10) ,1 

10 XI (I) =KSI (I) * *2 + (1. ODO-MACH(I) **2) **2 
X2 (I) =X1 (I) ** (1 . 0D0/3 . 0D0) 

XDOT ( I ) =1 . 0D0/X2 (I) 

GOTO 30 

20 XI (I) =- (KSI (I) * *2 - (1. ODO-MACH(I) **2) **2) 

X2 (I) =-Xl (I) ** (1 . 0D0/3 . 0D0) 

XDOT ( I ) =1 . 0D0/X2 (I) 

3 0 RETURN 
END 

c *********************************************************** 
C SUBROUTINE NEWTONCOTES * 

SUBROUTINE NEWTONCOTES (KSI, XDOT , X , S AVEX , SAVED , H, N, INDEX , I ) 
IMPLICIT REAL*8 (A-H,0-Z) 

IMPLICIT REAL* 8 (M) 

IMPLICIT REAL* 8 (K) 

DIMENSION X (N) ,XDOT(N) ,SAVEX(N) , SAVED (N) ,KSI (N) 

GO TO (1,2, 3, 4, 5) , index 

1 SAVEX ( I ) =X ( I ) 

SAVED ( I ) =XDOT ( I ) 

KSI (I) =KSI (I) +0 . 25D0*H 
RETURN 

2 SAVED ( I ) =7 . 0D0 * SAVED ( I ) +3 2 . ODO*XDOT ( I ) 

KSI (I) =KSI (I) +0 .25D0*H 

RETURN 

3 SAVED (I) = SAVED ( I ) +12 . ODO*XDOT ( I ) 

KSI (I) =KSI (I) +0 .25D0*H 

RETURN 

4 SAVED ( I ) = SAVED ( I) +32 . ODO*XDOT ( I) 

KSI (I) =KSI (I) +0 . 25D0*H 

RETURN 

5 X(I) =SAVEX ( I ) + (H/90.0D0) * (SAVED (I) +7 . ODO*XDOT (I) ) 

RETURN 

END 
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★ ********★★*★*★***★*★***★*★********★**★★★**★■*•***★*★*■*•■*■****** 



THIS FORTRAN PROGRAM USES NUMERICAL INTEGRATION BASED ON 
THE NEWTON-COTES FORMULA TO SOLVE NUMERICALLY EQNS . 41 AND 43. 

THE OUTPUT IS IN THE FORM OF ^ vs. Y 

program ETTA 
PARAMETER (N=5) 

PARAMETER (K=2000) 

IMPLICIT REAL*8 (A-H,0-Z) 

IMPLICIT REAL*8 (M) 

INTEGER INDEX, I, J 

DIMENSION Y (N) , YDOT (N) ,SAVEY(N) , SAVED (N) , TF (N) , MACH (N) 
DIMENSION ETA (N) , ETA11 (K) , Yll (K) , Y1 (N) , Y2 (N) 

MACH ( 1 ) =0 . 8D0 
MACH (2) =0 . 9D0 
MACH (3) =1 . 0D0 
MACH (4) =1 . 1D0 
MACH ( 5 ) =1 . 2D0 
H=0 . 01D0 

C INTEGRAL FINAL LIMITS * 

DO 7 1=1,5 

IF (MACH (I) .EQ.10D0) THEN 
TF ( I ) =0 . 0D0 
GOTO 7 
END IF 

IF (MACH ( I) . GT . 1 . 0D0 ) THEN 

TF ( I) = ( - (1 . 0D0 - MACH ( I ) **2) ) ** (1 . 0D0/3 . 0D0) 

GOTO 7 
END IF 

TF (I) = (1. ODO-MACH(I) **2) ** ( 1 . 0D0/3 . 0D0 ) 

7 END DO 

c MAIN PROGRAM * 

DO 100 1=1,5 
J=1 

ETA ( I ) =14 . 0D0 

Y(I)=2.0D0/ (ETA(I) ** ( 1 . 0D0/2 . 0D0 ) ) 

10 IF (ETA(I) -H.GT.TF (I) ) THEN 
20 do 50 index=l,5 

CALL DERIVative ( ETA, N, I , Y1 , Y2 , YDOT, MACH) 

CALL NEWTONCOTES ( ETA , YDOT , Y , S AVE Y , SAVED , H , N, INDEX , I ) 
50 end do 

ETA11 (J) =ETA ( I ) 

Yll (J) =Y ( I ) 

J=J+1 
GOTO 10 
END IF 
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H=0 . 001D0 

IF (ETA(I) -H.GT.TF (I) ) THEN 
GOTO 20 
END IF 

80 DO 90 L= J- 1,1, -1 

WRITE (20+i, *) Yll (L) , ETA11 (L) 

90 END DO 
100 END DO 
STOP 
END 

c SUBROUTINE DERIVative * 

£*★********★**★★★★★***★★★*★★*★***★★★★**★★************★★***** 
SUBROUTINE DERIVative (ETA, N , I , Y1 , Y2 , YDOT , MACH) 

IMPLICIT REAL* 8 (A-H, O- Z) 

IMPLICIT REAL* 8 (M) 

DIMENSION YDOT (N) , Y1 (N) , Y2 (N) , ETA ( N ) , MACH (N) 

GO TO (20,20,10,10,10) , I 

10 Y1 (I) -ETA (I) * *3 + ( 1 . 0D0 -MACH (I) **2) 

Y2 (I) =Y1 (I) ** (1.0D0/2 .0D0) 

YDOT ( I ) =1 .0D0/Y2 (I) 

GOTO 30 

2 0 Y1 (I) -ETA (I) * *3 - ( 1 . 0D0 -MACH ( I ) **2) 

Y2 ( I ) =Y1 ( I ) * * (1.0D0/2.0D0) 

YDOT (I) =1 . 0D0/Y2 (I) 

3 0 RETURN 

END 

c SUBROUTINE NEWTONCOTES * 

SUBROUTINE NEWTONCOTES ( ETA , YDOT , Y , SAVEY , SAVED , H, N, INDEX , I ) 
IMPLICIT REAL* 8 (A-H,0-Z) 

IMPLICIT REAL* 8 (M) 

DIMENSION Y (N) , YDOT (N) , SAVEY (N) , SAVED (N) , ETA ( N ) 

GO TO (1,2, 3, 4, 5) , index 

1 SAVEY (I) =Y (I) 

SAVED (I) = YDOT (I) 

ETA ( I ) =ETA (I) -0 .25D0*H 
RETURN 

2 SAVED (I) =7 . 0D0*SAVED (I) +32 . 0D0*YDOT ( I) 

ETA (I) =ETA ( I ) -0 .25D0*H 

RETURN 

3 SAVED (I) = SAVED ( I ) +12 . 0D0*YDOT(I) 

ETA ( I ) =ETA ( I ) -0 .25D0*H 

RETURN 

4 SAVED (I) = SAVED ( I ) +32 . 0D0*YDOT ( I) 

ETA ( I ) =ETA ( I ) - 0 . 25D0*H 

RETURN 

5 Y ( I ) =SAVEY ( I ) + (H/90.0D0) * (SAVED (I) +7 . 0D0 *YDOT ( I ) ) 

RETURN 

END 
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APPENDIX B 
FORM OF a AND /3 



The expressions for a and /3 as introduced in 

a = ±q (1 -Ml) 2 



[REF. 2] are 

(117) 



P = +C 2 (1 -Ml) (118) 

for 0.8SM.S1.2 . 

Upon inserting eqns . 117 and 118 in eqns . 16 and 21 yields 



X = f — dZ (119) 

o ^/z 2 + (1 -M.V 



-- / 

-(i-Ad) 3 



dw 

s[vP + (1-Mi) 



( 120 ) 



These relations are numerically integrated and plotted in 
Figs. 16 and 17. We observe the same solutions for M^sl.O as 
expected, but distinguish a different set of solutions for the 
subsonic case M oo <1.0 . 

Applying the same procedure as in Chap. IV, we easily 
compute the expression for the boundary surfaces 



dY 

\j fl 3 + (i -id) 



l ) 2 I i 

3 ' mZ( y+i 



( 121 ) 



37 



The numerical integration of eqn. 121 is shown in Fig. 18 for 
various M^, (e.g., 0.8,0.9,0.97,1.0,1.1,1.2). The diverging 
contours of the subsonic solutions appear to give 
contradictory results once investigating the pressure 
coefficient behavior 



C p = 



-2 

Af- ( Y + l) 



n 



2 + (l-M 2 ) 2 + (1-Mi) 



( 122 ) 



along the boundary surface. This implies as shown in Fig. 19 
that the boundary surface within the portion of C p >1.0 is not 
a realistic solution. 
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APPENDIX C 
FIGURES 



KSI(X) for M =0.8- 1.2 




Figure 1. Numerical integration of eqns. 40 and 42. 



ETTA(Y) for M =0.8- 1 .2 




Figure 2. Numerical integration of eqns. 41 and 43. 



E(Y) vs. Y for M =0.8- 1.2 




Figure 3. Numerical solution of E(Y). 



YO vs. M„ 




Figure 4. Equation 71. 



K(X) vs. X for M =0.8- 1.2 




Figure 5. Numerical integration of eqn. 72. 
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Numerical integration of eqn. 76. 
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Figure 7. Subsonic boundary surfaces. 
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Figure 8. Numerical solution of eqn. 77. 
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Figure 9. C p and local Mach for M CD =1.0 . 




Figure 10. C p and local Mach for M 00 =l.l . 
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Figure 11. C and local Mach for 
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Figure 12. C and local Mach for 
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Figure 13 . Velocity vectores for 




Figure 14. Mach lines for M tt =1.0 
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Figure 15. Upper transonic surface for 
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Figure 15b. Upper supersonic surface for 
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Figure 16. Numerical integration of eqn. 119. 



ETTA(Y) FOR MACH 0.8 - 1 .2 




Figure 17. Numerical integration of eqn. 120. 
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Figure 19. C p for M^O.8 . 
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